In [1]:
%matplotlib inline
import numpy as np
from datetime import datetime, timedelta
In [2]:
from siphon.catalog import TDSCatalog
# copied from the browser url box
metar_cat_url = 'http://thredds.ucar.edu/thredds/catalog/nws/metar/ncdecoded/catalog.xml?dataset=nws/metar/ncdecoded/Metar_Station_Data_fc.cdmr'
# parse the xml
metar_cat = TDSCatalog(metar_cat_url)
# what datasets are here? only one "dataset" in this catalog
dataset = list(metar_cat.datasets.values())[0]
print(dataset.name)
In [3]:
print(list(dataset.access_urls))
In [4]:
ncss_url = dataset.access_urls["NetcdfSubset"]
In [5]:
from siphon.ncss import NCSS
ncss = NCSS(ncss_url)
In [6]:
ncss.variables
Out[6]:
In [7]:
# get current date and time
now = datetime.utcnow() - timedelta(days=5)
now = datetime(now.year, now.month, now.day, now.hour)
# define the time range we are interested in
start_time = now - timedelta(days=1)
end_time = now
# build the query
query = ncss.query()
query.lonlat_point(-104.66, 39.85)
query.time_range(start_time, end_time)
query.variables('inches_ALTIM', 'air_temperature', 'dew_point_temperature',
'wind_from_direction', 'wind_speed')
query.accept('netcdf')
# what does the request url look like?
print(query)
In [8]:
data = ncss.get_data(query)
In [9]:
print(list(data.variables))
In [10]:
station_id = data['station_id'][:].tostring()
print(station_id)
In [11]:
station_id = station_id.decode("utf-8")
print(station_id)
In [12]:
time = [datetime.fromtimestamp(t) for t in data['time']]
In [13]:
from matplotlib import pyplot as plt
from matplotlib.dates import HourLocator,DateFormatter, AutoDateLocator
fig, ax1 = plt.subplots(1, 1)
ax1.plot(time, data['air_temperature'], '*')
locator = AutoDateLocator()
hoursFmt = DateFormatter('%H')
ax1.xaxis.set_major_locator(locator)
ax1.xaxis.set_major_formatter(hoursFmt)
ax1.autoscale_view()
ax1.set_title('Site: {} Date: {}'.format(station_id, time[0].strftime('%Y/%m/%d')))
ax1.set_xlabel('Hour of day')
#ax1.set_ylabel(t_air_label)
fig.autofmt_xdate()
plt.show()
So we have dewpoint, can we do a time series of mixing ratio?
We can use MetPy's unit support and calculations for this.
In [14]:
import metpy.calc
from metpy.units import units
In [15]:
print(data.variables['dew_point_temperature'].units)
print(data.variables['inches_ALTIM'].units)
In [16]:
dewp = data.variables['dew_point_temperature'][:] * units('degC')
slp = data.variables['inches_ALTIM'][:] * units('inHg')
Look at the docs for metpy.calc.mixing_ratio and metpy.calc.saturation_vapor_pressure.
So to get mixing ratio, we need the ambient partial pressure--we can get that by passing dewpoint temperature to metpy.calc.saturation_vapor_pressure().
In [17]:
e = metpy.calc.saturation_vapor_pressure(dewp)
e
Out[17]:
In [18]:
slp
Out[18]:
So we can pass these to mixing ratio and get...
In [19]:
m = metpy.calc.mixing_ratio(e, slp)
m
Out[19]:
Hmmm....need to get this to reduce
In [20]:
# Could also to m.ito('dimensionless')
m.ito_base_units()
m
Out[20]:
In [21]:
from matplotlib import pyplot as plt
from matplotlib.dates import HourLocator,DateFormatter, AutoDateLocator
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(time, data['air_temperature'], 'r*')
ax1.plot(time, data['dew_point_temperature'], 'go')
ax1.grid(True)
ax2.plot(time, m, 'x')
ax2.grid(True)
locator = AutoDateLocator()
hoursFmt = DateFormatter('%H')
ax1.xaxis.set_major_locator(locator)
ax1.xaxis.set_major_formatter(hoursFmt)
ax1.autoscale_view()
ax1.set_title('Site: {} Date: {}'.format(station_id, time[0].strftime('%Y/%m/%d')))
ax1.set_xlabel('Hour of day')
fig.autofmt_xdate()
plt.show()
In [ ]:
In [22]:
bb = {'north' : 45,
'south' : 35,
'east' : -100,
'west' : -110}
query = ncss.query()
query.lonlat_box(north=bb['north'], south=bb['south'], east=bb['east'], west=bb['west'])
query.time(start_time)
query.variables('air_temperature', 'dew_point_temperature', 'inches_ALTIM',
'wind_speed', 'wind_from_direction', 'cloud_area_fraction')
query.accept('csv')
Out[22]:
In [23]:
from metpy.calc import get_wind_components
data = ncss.get_data(query)
# Access is just like netcdf4-python
lats = data['latitude'][:]
lons = data['longitude'][:]
tair = data['air_temperature'][:]
dewp = data['dew_point_temperature'][:]
slp = (data['inches_ALTIM'][:] * units('inHg')).to('mbar')
# Convert wind to components
u, v = get_wind_components(data['wind_speed'], np.deg2rad(data['wind_from_direction']))
# Need to handle missing (NaN) and convert to proper code
cloud_cover = 8 * data['cloud_area_fraction']
cloud_cover[np.isnan(cloud_cover)] = 9
cloud_cover = cloud_cover.astype(np.int)
# For some reason these come back as bytes instead of strings
stid = [s.decode() for s in data['station']]
Pull out the code to create the basic plot and add mapping features so we aren't repeating this everywhere.
In [24]:
import cartopy
def default_map(bb):
fig = plt.figure(figsize=(24, 12))
proj = cartopy.crs.Stereographic(central_longitude=-95, central_latitude=35)
ax = fig.add_subplot(1, 1, 1, projection=proj)
# Add map features
ax.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural',
name='admin_1_states_provinces_lakes',
scale='50m',
facecolor='none'))
ax.add_feature(cartopy.feature.BORDERS)
ax.coastlines()
ax.gridlines()
# Set extent to match requested bounding box
ax.set_extent([bb['west'], bb['east'], bb['south'], bb['north']])
return ax
In [25]:
from metpy.plots import StationPlot
from metpy.plots.wx_symbols import sky_cover
ax = default_map(bb)
# Start the station plot by specifying the axes to draw on, as well as the
# lon/lat of the stations (with transform). We also set the fontsize.
stationplot = StationPlot(ax, lons, lats, transform=cartopy.crs.PlateCarree(),
fontsize=12)
# Plot the temperature and dew point to the upper and lower left, respectively, of
# the center point. Each one uses a different color.
stationplot.plot_parameter('NW', tair, color='red')
stationplot.plot_parameter('SW', dewp, color='darkgreen')
# Add wind barbs
stationplot.plot_barb(u, v)
In addition to plotting values, StationPlot has support for plotting text strings, symbols, and plotting values using custom formatting.
Plotting symbols involves mapping integer values to various custom font glyphs in our custom weather symbols font. MetPy provides mappings for converting WMO codes to their appropriate symbol. The sky_cover function below is one such mapping.
Below we also use a custom formatter to take the sea level pressure values and plot them in the prototypical 3 digit form.
In [26]:
from metpy.plots.wx_symbols import sky_cover
ax = default_map(bb)
# Same as before
stationplot = StationPlot(ax, lons, lats, transform=cartopy.crs.PlateCarree(),
fontsize=12)
stationplot.plot_parameter('NW', tair, color='red')
stationplot.plot_parameter('SW', dewp, color='darkgreen')
stationplot.plot_barb(u, v)
# Plot the sky cover symbols in the center. We give it the integer code values that
# should be plotted, as well as a mapping class that can convert the integer values
# to the appropriate font glyph.
stationplot.plot_symbol('C', cloud_cover, sky_cover)
# Plot station id -- using an offset pair instead of a string location
stationplot.plot_text((2, 0), stid)
# A more complex example uses a custom formatter to control how the sea-level pressure
# values are plotted. This uses the standard trailing 3-digits of the pressure value
# in tenths of millibars.
stationplot.plot_parameter('NE', slp,
formatter=lambda v: format(10 * v.magnitude, '.0f')[-3:])
Out[26]:
Station plots can also be created using layouts, which encapsulate the formatting elements for various data types (based on string name). Formatting elements include:
In [27]:
from metpy.plots import simple_layout
simple_layout.names()
Out[27]:
Using the layout takes data arrays stored in a dictionary like object, where the keys represent data fields. (In an ideal world, you could just pass a netCDF4 Dataset instance with CF-compliant names...)
The layout will ignore any fields that aren't present. When given data, it will ignore any names (keys) that are not specified in the layout. It is designed to be as forgiving as possible.
In [28]:
sfc_data = {'air_temperature': tair, 'dew_point_temperature': dewp, 'eastward_wind': u,
'northward_wind': v, 'cloud_coverage': cloud_cover,
'air_pressure_at_sea_level': slp}
In [29]:
ax = default_map(bb)
stationplot = StationPlot(ax, lons, lats, fontsize=12,
transform=cartopy.crs.PlateCarree())
simple_layout.plot(stationplot, sfc_data)
You can also create your own layout:
In [30]:
from metpy.plots import StationPlotLayout
layout = StationPlotLayout()
layout.add_barb('u', 'v', 'knots')
layout.add_symbol('C', 'sky', sky_cover)
# These are wider fields, so we'll put them out wider
layout.add_text((2, 0), 'station', color='blue')
layout.add_value((-2, 0), 'temp', fmt='0.1f', units='degF', color='red')
So we'll put data into a dictionary to display...
In [31]:
sfc_data = {'temp': tair * units('degC'), 'sky':cloud_cover, 'u': u, 'v': v,
'station': stid}
and create the plot.
In [32]:
ax = default_map(bb)
stationplot = StationPlot(ax, lons, lats, fontsize=12,
transform=cartopy.crs.PlateCarree())
layout.plot(stationplot, sfc_data)
MetPy has many more symbols for current weather and cloud types; we just lack the data to readily show them off at present. Exhaustive list:
These all assume WMO code values for these items, and they all can be passed as a mapper for a symbol field in the station plot code (either plot_symbol for StationPlot or add_symbol for StationPlotLayout)
Below we show code to display all of the weather symbols:
In [33]:
# Get the mapper that converts a code for current weather (from WMO) to the appropriate
# unicode code point
from metpy.plots.wx_symbols import current_weather, wx_symbol_font
# Need to make a new copy of the font so we can make it bigger
big_font = wx_symbol_font.copy()
big_font.set_size(36)
# Create a plot to loop over all of the codes and display it
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
for i in range(100):
ax.text(i % 10, i // 10, current_weather(i), fontproperties=big_font)
ax.set_xlim(0, 10)
ax.set_ylim(0, 10);
Or sky cover symbols:
In [34]:
fig, ax = plt.subplots(1, 1, figsize=(8, 2))
for i in range(10):
ax.text(i, 0, sky_cover(i), fontproperties=big_font)
ax.set_xlim(0, 10)
ax.set_ylim(0, 1);
Time to play around with the station plot. Take some of the data we grabbed (or change the query and download data for a different area) and pass it to the station plotting routines. You can use either the direct StationPlot class or create your own layout. If you want to play with the various symbols, you can randomly generate integer codes using numpy.random.randint(low, high)
In [ ]: